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The Tokai-to-Kamioka (T2K) experiment studies neutrino oscillations using an off-axis muon 
neutrino beam with a peak energy of about 0.6 GeV that originates at the J-PARC accelerator 
facility. Interactions of the neutrinos are observed at near detectors placed at 280 m from the 
production target and at the far detector - Super-Kamiokande (SK) - located 295 km away. The 
flux prediction is an essential part of the successful prediction of neutrino interaction rates at the T2K 
detectors and is an important input to T2K neutrino oscillation and cross section measurements. 
A FLUKA and GEANT3 based simulation models the physical processes involved in the neutrino 
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production, from the interaction of primary beam protons in the T2K target, to the decay of hadrons 
and muons that produce neutrinos. The simulation uses proton beam monitor measurements as 
inputs. The modeling of hadronic interactions is re-weighted using thin target hadron production 
data, including recent charged pion and kaon measurements from the NA61/SHINE experiment. 
For the first T2K analyses the uncertainties on the flux prediction are evaluated to be below 15% 
near the flux peak. The uncertainty on the ratio of the flux predictions at the far and near detectors 
is less than 2% near the flux peak. 

PACS numbers: 24.10.Lx,14.60.Lm 



I. INTRODUCTION 

Predicting the neutrino flux and energy spectrum is an 
important component of analyses in accelerator neutrino 
experiments [1-4]. However, it is difficult to simulate 
the flux precisely due to uncertainties in the underly- 
ing physical processes, particularly hadron production 
in proton- nucleus interactions. To reduce flux-related 
uncertainties, neutrino oscillation experiments are some- 
times conducted by comparing measurements between a 
near detector site and a far detector site, allowing for 
cancellation of correlated uncertainties. Therefore, it is 
important to correctly predict the relationship between 
the fluxes at the two detector sites, described below as 
the far-to-near ratio. 

T2K (Tokai-to-Kamioka) [5] [6] is a long-baseline neu- 
trino oscillation experiment that uses an intense muon 
neutrino beam to measure the mixing angle #13 via the 
v e appearance [7] and the mixing angle #23 and mass dif- 
ference Am| 2 via the disappearance [8]. The muon 
neutrino beam is produced as the decay products of pi- 
ons and kaons generated by the interaction of the 30 GeV 
proton beam from Japan Proton Accelerator Research 
Complex (J- PARC) with a graphite target. The prop- 
erties of the generated neutrinos are measured at near 
detectors placed 280 m from the target and at the far 
detector, Super-Kamiokande (SK) [9], which is located 
295 km away. The effect of oscillation is expected to be 
negligible at the near detectors and significant at SK. 

The T2K experiment employs the off-axis method [10] 
to generate a narrow-band neutrino beam and this is the 
first time this technique has been used in a search for neu- 
trino oscillations. The method utilizes the fact that the 
energy of a neutrino emitted in the two-body pion (kaon) 
decay, the dominant mode for the neutrino production, 
at an angle relative to the parent meson direction is only 
weakly dependent on the momentum of the parent. The 
parent 7r + (~)'s are focused parallel to the proton beam 
axis to produce the (anti-)neutrino beam. By position- 
ing a detector at an angle relative to the focusing axis, 
one will, therefore, see neutrinos with a narrow spread 
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in energy. The peak energy of the neutrino beam can be 
varied by changing the off-axis angle as illustrated in the 
lower panel of Fig. 1. In the case of T2K, the off-axis 
angle is set at 2.5° so that the neutrino beam at SK has 
a peak energy at about 0.6 GeV, near the expected first 
oscillation maximum (Fig. 1). This maximizes the effect 
of the neutrino oscillations at 295 km as well as reduces 
background events. Since the energy spectrum changes 
depending on the off-axis angle, the neutrino beam di- 
rection has to be precisely monitored. 
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E v (GeV) 

FIG. 1: Muon neutrino survival probability at 295 km 
and neutrino fluxes for different off-axis angles. 

To determine the oscillation parameters, the expected 
observables at the far detector are predicted based on 
the flux prediction and the neutrino-nucleus interaction 
model. To reduce the uncertainty of the prediction, they 
are modified based on the near detector measurements. 
For example, the absolute normalization uncertainty is 
efficiently canceled by normalizing with the event rate at 
the near detector. Then, it is important to reduce the 
uncertainty on the relation between the flux at the near 
detector and that at the far detector. 

The physics goals of T2K are to be sensitive to the val- 
ues of sin 2 26*13 down to 0.006 and to measure the neu- 
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trino oscillation parameters with precision of 5(Am| 2 ) ~ 
10" 4 eV 2 and 5(sin 2 20 2 3) ~ 0.01. To achieve these, the 
near-to-far extrapolation of the flux, i.e., the far-to-near 
flux ratio as a function of energy has to be known to bet- 
ter than 3%. In addition to this requirement, it is also 
desirable to reduce the absolute flux uncertainty to study 
the neutrino-nucleus interactions at the near detector. 

For this purpose, the fluxes are calculated and the 
uncertainties are estimated based on hadron production 
measurements including those by the NA61/SHINE ex- 
periment [11] [12] and in situ measurements of the pri- 
mary proton beam properties and the neutrino beam di- 
rection. 

In this paper, we describe a Monte Carlo based neu- 
trino flux prediction as a function of neutrino energy at 
near and far detectors in T2K and the methods to esti- 
mate the flux prediction uncertainties. The neutrino flux 
treated here is the flux for the 'neutrino' running mode, in 
which positive pions are focused. Section II describes the 
neutrino beamline, while Sec. Ill summarizes the beam 
operation history. Section IV describes a method of neu- 
trino flux prediction based on a data-driven simulation. 
Section V explains uncertainties on the flux prediction. 
A comparison between the measured and predicted flux 
is discussed in Sec. VI. 




(3) Final focusing section I 

(4) Target station i 

(5) Decay volume J , 

(6) Beam dump j 

FIG. 2: An overview of the T2K neutrino beamline. 



five current transformers (CTs), 21 electrostatic moni- 
tors (ESMs), 19 segmented secondary emission monitors 
(SSEMs), respectively. The monitor locations in FF sec- 
tion are shown in Fig. 3. 



II. T2K NEUTRINO BEAMLINE 

The J-PARC Main Ring (MR) accelerates a 30 GeV 
proton beam every 2 to 3 seconds. For each acceleration 
cycle, the beam is fast-extracted to the T2K neutrino 
beamline as a 'spill'. One spill contains eight bunches in 
about 5 /is. 

The neutrino beamline is composed of two sections: 
the primary and secondary beamlines. In the primary 
beamline, the extracted proton beam is transported to 
point in the direction of the secondary beamline, and fo- 
cused to have the desired profile at the target. In the 
secondary beamline, the proton beam impinges on a tar- 
get to produce secondary pions and other hadrons, which 
are focused by magnetic horns and decay into neutrinos. 
An overview of the neutrino beamline is shown in Fig. 2. 
More details of the beamline are described in [6]. 



A. Primary beamline 

The primary beamline consists of the preparation sec- 
tion (54 m long), arc section (147 m) and final focus- 
ing section (37 m). In the final focusing (FF) section, 
ten normal conducting magnets (four steering, two dipole 
and four quadrupole magnets) guide and focus the beam 
onto the target, while directing the beam downward by 
3.64 degrees with respect to the horizontal. 

The intensity, position and profile of the proton beam 
in the primary sections are precisely monitored by 



1. Proton Beam Monitor 

The beam intensity is measured with five CTs. Each 
CT is a 50-turn toroidal coil around a cylindrical ferro- 
magnetic core. The uncertainty on the beam intensity 
is 2%, which originates from the calibration accuracy 
(1.7%), the effect of secondary electrons produced at the 
SSEM foils (<0.7%), the long term stability of the in- 
dividual CT monitors relative to each other and the CT 
monitor measurement from the main ring (0.5%). For the 
flux prediction, the intensity measured by CT5, located 
most downstream, is used. 

The ESMs have four segmented cylindrical electrodes 
surrounding the proton beam orbit. By measuring 
the top-bottom and left-right asymmetry of the beam- 
induced current on the electrodes, they monitor the pro- 
ton beam center position nondestructively (without di- 
rectly interacting with the beam). The measurement 
precision of the projected beam position is better than 
450 /zm. 

The SSEMs have two 5 /im thick sets of titanium foil 
strips oriented horizontally and vertically in the plane 
perpendicular to the beam axis, and a high voltage an- 
ode foil between them. They measure the horizontal and 
vertical projections of the proton beam profile. The sys- 
tematic uncertainty of the beam width measurement is 
200 /xm. The uncertainty of the beam center position 
measurement is dominated by the monitor alignment un- 
certainty discussed in Section II C. Since each SSEM 
causes a beam loss (0.005% loss), they are inserted into 
the beam orbit only during the beam tuning, and re- 
moved from the beam orbit during the continuous beam 
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FIG. 3: Location of the primary beamline monitors in the final focusing section. 



operation except for the most downstream SSEM. 

An optical transition radiation (OTR) monitor posi- 
tioned 30 cm upstream of the target measures the two 
dimensional profiles of the beam by imaging transition 
radiation produced when the beam crosses a 50 /im thick 
titanium alloy foil. The details of the monitor have been 
described elsewhere [13]. 

Using the ESMs, SSEMs and OTR measurements, the 
beam position at the upstream side of the baffle (shown in 
Fig. 4) is reconstructed with accuracy better than 0.7 mm 
as described in Sec. Ill A. 



B. Secondary beamline 
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Pions and kaons are produced by the interaction of 
protons with a graphite target. They decay in-flight in- 
side a single volume of ^1500 m 3 filled with helium gas. 
The helium vessel is connected with the primary beam- 
line using a titanium-alloy beam window that separates 
the vacuum in primary beamline and helium gas volume 
in the secondary beamline. 

The secondary beamline consists of three sections: the 
target station, decay volume and beam dump (Fig. 4). 
The helium vessel in the target station is 15 m long, 4 m 
wide and 11 m high. The decay volume is a 96 m long 
steel tunnel. The cross section is 1.4 m wide and 1.7 m 
high at the upstream end, and 3.0 m wide and 5.0 m high 
at the downstream end. The beam dump sits at the end 
of the decay volume. The distance between the center of 
the target and the upstream surface of the beam dump 
is 109 m. 

The target station contains a baffle, the OTR moni- 
tor, the target and three magnetic horns. The baffle is 
a collimator to protect the horns. The 250 kA current 
pulses magnetize the three horns to focus the secondary 
7r + 's in 'neutrino' running mode. The 7r _ 's are focused in 
'anti-neutrino' running mode, where the polarity of the 
horn current is inverted. The produced pions then decay 
in the decay volume mainly into muons and muon neu- 
trinos. All the remnants of the decayed pions and other 
hadrons are stopped by the beam dump. The neutrinos 
pass through the beam dump and are used for physics 
experiments. The muons above 5 GeV that also pass 



FIG. 4: Side view of the secondary beamline. 



through the beam dump are detected by a muon mon- 
itor (MUMON) that monitors the beam direction and 
intensity. 



1. Target and Horns 

The target core is a 1.9 interaction length (91.4 cm 
long), 2.6 cm diameter graphite rod with a density of 1.8 
g/cm 3 . The core and a surrounding 2 mm thick graphite 
tube are sealed inside a 0.3 mm thick titanium case. The 
target assembly is cantilevered inside the bore of the first 
horn inner conductor. 

T2K uses three magnetic horns. Each horn consists 
of two coaxial (inner and outer) conductors which en- 
compass a closed volume [14, 15]. A toroidal magnetic 
field is generated in that volume. The field varies as 
1/r, where r is the distance from the horn axis. The 
first horn (Horn 1) collects the pions that are generated 
at the target installed in its inner conductor. The sec- 
ond (Horn 2) and third (Horn 3) horns focus the pions. 
When the horns are operating with a current of 250 kA, 
the maximum field is 1.7 T and the neutrino flux at SK is 
increased by a factor of ~17 at the spectrum peak energy 
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TABLE I: Dimensions of the T2K horns 
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FIG. 5: The predicted flux of at the SK far detector 
for operation at different horn currents. The flux 
histogram (top) ranges from to 3 GeV, while the 
ratios (bottom) range from to 10 GeV. 
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FIG. 6: Cross section view of horns. 



(~0.6 GeV), as illustrated in Fig. 5. 

A schematic view of the horns is shown in Fig. 6. The 
horn conductor is made of an aluminum alloy. Their 
dimensions are summarized in Table I. The thickness 
of the inner conductors is 3 mm. They are optimized 
to maximize the neutrino flux; the inside diameter is as 
small as possible to achieve the maximum magnetic field, 
and the conductor is as thin as possible to minimize pion 
absorption while still being tolerant of the Lorentz force 
from the 250 kA current and the thermal shock from the 
beam [16]. 

The electrical currents of the magnetic horns are mon- 



Hornl Horn2 Horn3 



inner conductor inside diameter (mm) 54 80 140 
outer diameter (mm) 400 1000 1400 

length (m) 1.5 2 2.5 



TABLE II: Uncertainties on the absolute horn current 
measurement. In the total error calculation, full width 
(FW) errors are scaled by l/\/l2 to estimate la- 
uncertainty. 



uncertainty 


coil calibration 


±1% (FW) 


coil setting 


±1% (FW) 


electronics calibration 


< 1% 


monitor stability 


2% (FW) 


total 


1.3% 



itored by Rogowski coils whose signal are digitized by 
65 MHz FADCs. Table II shows the summary of the 
horn current uncertainties. The Rogowski coils were cal- 
ibrated by the production company with ±1% precision. 
The shape of the "loop" of the Rogowski coil may cause 
a 1% change of gain. 

FADCs and related electronics are calibrated with bet- 
ter than 1% precision. 

Each horn has several instrumentation ports at vari- 
ous positions along the horn axis which permit measure- 
ments of the magnetic field between the inner and outer 
conductors. Multiple magnetic field measurements have 
been made on the horns to validate the nominal 1/r field 
and to check for the presence of magnetic field asymme- 
tries. The magnetic fields generated by Horns 2 and 3 
were measured using an integrated 3-axis Hall probe in- 
serted between the inner and outer conductors via the 
horns' instrumentation ports. The results are summa- 
rized in Table III. The measured field agrees with the 
expected nominal field within 2%. 

Measurements of the magnetic field were also taken 
on a spare copy of the first horn, identical in design 
to the one currently in use in the T2K beamline. As 
with Horns 2 and 3, field measurements were taken via 
the instrumentation ports using a 3-axis Hall probe. A 
comparison of the expected field to the data taken at 
the right upstream port is shown in Fig. 7. The results 
agree well with the expected nominal field. Additional 
measurements were taken along the horn's axis inside 
of the inner conductor. The purpose of these measure- 
ments was to detect possible magnetic field asymmetries 
caused by path length differences between the upper and 
lower striplines supplying current to the horn. While 
no field asymmetry due to path length differences was 
observed, an on-axis magnetic field with an anomalous 
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FIG. 7: Measurements of the magnetic field magnitude 
taken at the right upstream port of Horn 1. The curve 
shows the expected field strength, including a small 
correction to account for fringe effects near the 
instrumentation port at large radii. 

TABLE III: Magnetic field deviations from expected 
values at all instrumentation ports. Blanks in the table 
are a result of each horn having a different configuration 
of instrumentation port locations. 



Top 


Bottom 


Left 


Right 


Horn 1 Upstream 

Downstream 




0.94% 


0.5% 
1.0% 


Horn 2 Upstream 0.7% 
Midstream 0.7% 


0.1% 
0.6% 


1.3% 


0.7% 


Horn 3 Upstream 1.2% 
Downstream 0.7% 




1.2% 
0.2% 


1.0% 
0.5% 



time-dependence was detected. While the magnitude of 
the azimuthal fields is always proportional to the current, 
the anomalous on-axis field is characterized by a differ- 
ence of 0.7 ms between maximum field and maximum 
current, as shown in Fig. 8. This field has a direction 
perpendicular to the beam axis, and reaches a maximum 
magnitude of 0.065 T near the center of the horn. The 
cause of this anomalous on-axis field is not yet known. 
Therefore, the effect of this field is estimated and added 
to the flux uncertainty. The magnitude of this field is 
3.7% of the magnitude of the field just outside the inner 
conductor, and conservative estimates from Monte Carlo 
simulations show the effect on neutrino flux is small, as 
discussed in VD. 



2. Secondary beam monitoring 

The neutrino beam intensity and direction can be mon- 
itored on a bunch-by-bunch basis by measuring the pro- 
file of muons which are produced along with neutrinos 
from the pion two-body decay. The neutrino beam direc- 
tion is measured as the direction from the target to the 



0.07 r 
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0.05- 
0.04^ 
0.03^ 
0.02^ 
0.01^ 



Measured field strength 
Measured current (not to scale) 



0123456789 10 

Time (\ls) 

FIG. 8: A sample of data from the Hall probe showing 

the field strength in the x direction in the beam 
coordinates. This data was taken 100 cm along the axis 

of the horn. The Rogowski coil output, showing the 
current time dependence, is drawn with unfilled markers 
and is not to scale. The peaks arc offset by 
approximately 0.7 ms. 



center of the muon profile. The muon monitor is located 
just behind the beam dump at a distance of 118 m from 
the target, as shown in Fig. 4. It consists of two kinds 
of detector arrays: ionization chambers and silicon PIN 
photodiodes. Each array consists of 49 sensors at 25 cm 
intervals and covers a 150 x 150 cm 2 area. The precision 
on the center of the muon profile is 2.95 cm, which cor- 
responds to 0.25 mrad precision on the beam direction. 
The details of this monitor are described in [17]. 

The neutrino beam intensity and direction are mon- 
itored directly by measuring the profile of neutrinos at 
the INGRID detector [18], located 280 m away from the 
target. It consists of 16 identical neutrino detectors ar- 
ranged in horizontal and vertical arrays around the beam 
center. Each neutrino detector has a sandwich struc- 
ture of the iron target plates and scintillator trackers. 
The intensity and profile of the neutrino beam are re- 
constructed from the number of detected neutrino inter- 
actions in each module. At the present beam intensity 
(about 10 18 protons-on-target/day), the neutrino event 
rate is monitored daily with 4% precision. The neutrino 
beam center is measured monthly with accuracy better 
than 0.4 mrad. 

The ND280 detector measures the off-axis neutrino 
flux at a baseline of 280 m. At 280 m, ND280 ef- 
fectively sees a line source of neutrinos rather than a 
point source, therefore it covers a range of off-axis an- 
gles. The off-axis angle to ND280 from the target posi- 
tion is 2.04°. This angle was chosen to make the neu- 
trino spectrum at ND280 as similar as possible to the 
spectrum at SK. Consisting of electromagnetic calorime- 
ters, scintillating trackers and time projection chambers 
in a magnetic field, the ND280 detector can measure the 
spectrum of neutrino interactions for use in the extrapo- 
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lation of the flux prediction to SK. Independent neutrino 
cross section measurements can also be made at ND280, 
for which well-controlled absolute flux uncertainty is a 
key ingredient. The details of the ND280 detector arc 
given in [6, 19, 20]. 



C. Alignment of the beamline and actual neutrino 
beam direction 

The neutrino beam direction with respect to SK ( "off- 
axis angle" ) and the distance between the target and SK 
are obtained by GPS survey. The distance between the 
target and the center position of SK is 295,335.9±0.7 m. 
The beam axis is declined by 3.637°, while SK direction is 
1.260° downward and 0.795° to the north from the beam 
axis. The off-axis angle is adjusted to 2.50° to maximize 
the neutrino oscillation probability and measured to be 
2.504±0.004°. 

Based on the surveys, the primary beamline compo- 
nents, target, and horns were aligned in order to send the 
neutrino beam in the right direction. The muon moni- 
tor and the neutrino near detectors were also aligned in 
order to monitor the neutrino beam direction. 

The directional accuracy of a long-baseline GPS survey 
is about 3xl0~ 6 rad. The accuracy of a short distance 
survey (e.g. the direction to the near detectors) is about 
7xl0~ 5 rad. 

The alignment of the components in the primary beam- 
line was carried out with a laser tracker which has a 
spatial resolution of 50 /zm for distances shorter than 
20 m. The proton monitors were aligned to better than 
0.4 mm. The OTR monitor, in the secondary beamline, 
was aligned with a precision of 0.3 mm relative to the 
axis of the first horn. The relative alignment between the 
OTR and upstream proton monitors is known to within 
1.0 mm. 

The target was surveyed relative to the horn axis after 
installation. A difference of 1.2 mm (0.3 mm) in hor- 
izontal (vertical) direction between the target and the 
horn axis was measured at the downstream end, while 
the alignment of the upstream end was found to be cen- 
tered on the horn axis to within 0.1 mm. 

The observed displacement at the downstream end of 
the target translates into 1.3 mrad (0.3 mrad) rotation 
in the horizontal (vertical) plane of the downstream end 
relative to the target head. The effect of this rotation on 
the predicted neutrino flux is included as a systematic 
uncertainty (see Section VC). 

The position of each horn was surveyed after the in- 
stallation. In the directions horizontally transverse and 
parallel (x and z) to the proton beam axis, the horns were 
aligned relative to the beamline survey markers inside the 
helium vessel. The alignment accuracy in these directions 
are 0.3 mm and 1.0 mm for x and z, respectively. The 
vertical position, y, of the center of each horn was aligned 
relative to the survey markers on one of the magnets in 
the final section of the primary beamline. The alignment 



precision in this direction is dominated by an overall un- 
certainty of 1.0 mm in the vertical alignment between 
the primary and secondary beamlines. The precision of 
the angular alignment of each horn is about 0.2 mrad, 
which is based on the survey of the alignment markers 
at the upstream and downstream end of each horn. The 
movement of the horn conductors by the current pulse 
was measured to be less than 0.1 mm. 

After the 2011 Tohoku earthquake, movements of the 
GPS survey points, of the primary beamline tunnel and 
of the beamline components were observed. The baseline 
to SK was increased by 0.9 m, while the beam angle was 
rotated by 3 x 10 -5 rad. Both of these shifts have a small 
effect on the physics performance of the experiment. The 
upstream end of the primary beamline tunnel was ob- 
served to have shifted by 12 mm horizontally and 5 mm 
vertically relative to its downstream end, which is fixed to 
the target station and secondary beamline. The beam- 
line magnets and monitors were realigned to the same 
alignment accuracy with the fixed point of reference at 
the most downstream end of the primary beamline. 

The horns were observed to have shifted by 3~9 mm 
according to the survey of alignment markers at the top of 
the horn support modules. The horns were realigned us- 
ing the survey markers on the support modules, and the 
alignment was confirmed by lowering a rigid frame with 
a camera and alignment laser into the helium vessel and 
checking the position of survey marks on the horns. The 
horns were found to be at the expected position, within 
the 1 mm accuracy of the survey method. The alignment 
of the OTR monitor could not be directly checked since 
the first horn was not removed from the helium vessel. 
In situ scans of the beam across the target, after realign- 
ment of the primary beamline monitors, have shown that 
the measured beam position by the OTR monitor is con- 
sistent with the beam position extrapolated from the pri- 
mary beamline SSEM monitors, as shown in Fig. 9. The 
MUMON was surveyed after the earthquake and its po- 
sition relative to the target station shifted by less than 
1 mm. 



III. BEAM OPERATION HISTORY 

The proton beam profile, neutrino beam direction and 
horn current are measured during the beam operation 
and measurement results are reflected in the neutrino 
flux prediction and estimation of uncertainty. The flux 
prediction and uncertainty presented in this paper are 
based on three physics runs: Run 1 (January - June 
2010), Run 2 (November 2010 - March 2011), Run 3 
(March - June 2012). The Run 3 period is divided into 
three sub periods - Run 3a (March 2012), Run 3b (March 
2012) and Run 3c (April - June 2012) - according to the 
horn current settings (with a kA setting in Run 3a 
and a 205 kA setting in Run 3b instead of the nominal 
250 kA). The polarity of the horn current was always 
set to that for the 'neutrino' running mode. The Run 
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FIG. 9: The correlations between the beam position 
measurements in x (top) and y (bottom) by the OTR 
monitor and the SSEM monitors extrapolated to the 
OTR position. The intercept and slope are from a 
linear fit (red line) to the measurements. 



3a data is not used for the oscillation analysis because 
the data in this period is small (0.3% of the total) and 
the horn current was set to kA. However, it is used for 
studies of kA operation. Figure 10 shows the plot of the 
accumulated POT and protons per pulse for good quality 
beam data over time. The total accumulated number of 
protons (protons on the target, POT) in all run periods is 
3.04x 10 20 POT, corresponding to 4% of T2K's exposure 
goal. The maximum beam power reached so far is about 
200 kW. 

We select only good quality beam data for physics anal- 
ysis using the following conditions. 

• Each hardware component works normally. 

• The deviation of all Horns currents from the mean 
is within ± 5 kA. 

• The deviation of the beam angle measured by MU- 
MON from the mean is within 1 mrad. 



• The deviation of the total muon yield measured by 
MUMON from the mean is within ± 5 %. 

The beam data from the beam monitors are checked on a 
spill-by-spill basis during beam operation to ensure data 
from good quality beam is used. For example, Fig. 11 
shows the history of the muon profile center measured 
at MUMON. In all run periods, this profile center is sta- 
ble within 1 mrad from the beam axis (1 mrad stability 
is required to achieve the physics goal of T2K). During 
Run 3b period, the MUMON center deviated from the 
beam-axis in both the X and Y directions. The possible 
reason for this deviation is misalignment of the horns. 
The beam may be focused in the deviated direction if 
there is a horn misalignment. This deviation can change 
depending on the horn current. As described later, the 
direction of the neutrino beam in this period had been 
also measured at INGRID and it was confirmed to be 
within 1 mrad (see Table VIII). After the good quality 
cut, the fraction of beam data is 99.8%. 

A. Proton beam properties 

The center position and angle of the proton beam at 
the upstream surface of the baffle are reconstructed by 
extrapolating the center positions measured with ESM20, 
SSEM19 and OTR for the vertical and ESM19, ESM20, 
SSEM19 and OTR for the horizontal direction for each 
spill. 

Each time the beam conditions change, all of the 
SSEMs are inserted into the beamline and beam pro- 
files are measured for 100 spills. The Twiss parameters 
and emittance are obtained by fitting these profiles along 
with that from the OTR. The beam width and diver- 
gence at the baffle are calculated from the Twiss param- 
eters and emittance. After 100 spills, all SSEMs except 
for SSEM19 are extracted from the beam orbit and the 
beam width and divergence are then obtained by scaling 
the emittance from the change of the profile measured at 
SSEM19 and OTR for each spill. 

Proton beam properties for each run period are ob- 
tained by reconstructing the beam profile at the baffle for 
each spill and summing the spill-by-spill profiles weighted 
by the number of protons in each spill. Tables IV and 
V summarize the measured center position, angle, width, 
emittance and Twiss a at the baffle for each run period. 

The following are the sources of uncertainty on the 
measurements of the average position and angle of the 
beam: 

• The alignment uncertainty of the proton beam 
monitors 

• The alignment uncertainty between the primary 
proton beamline and the target station (secondary 
beamline) 

• Systematic errors in the position measurements by 
the proton beam monitors 
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TABLE IV: Summary of measured proton beam 
parameters in the horizontal direction at the baffle for 
each run period : center position (X) and angle (Ox), 

Gaussian width (a), emittance (e), and Twiss (a). 



Period X (mm) 


6x (mrad) 


a (mm) 


e (ir mm mrad) 


a 


Runl 0.37 


-0.044 


4.27 


2.12 


0.60 


Run2 0.149 


-0.080 


4.04 


5.27 


0.16 


Run3b 0.087 


0.020 


4.13 


6.50 


0.16 


Run3c -0.001 


0.032 


4.03 


4.94 


0.33 



This is derived from the difference in the field inte- 
gral of dipole magnets between the estimate using 
the actual beam orbit and one obtained by mul- 
tiplying the measured field strength and the pole 
length. 

The resulting uncertainties on the emittance, Twiss a 
and beam width are summarized in Table VII. 

TABLE VII: Uncertainties for the emittance e, Twiss a 
and width a at baffle of the proton beam. 



TABLE V: Summary of measured proton beam 
parameters in the vertical direction at the baffle for 
each run period : center position (Y) and angle (0y), 
Gaussian width (a), omittance (e), and Twiss (a). 



ex f-Y ctx ay ox cry 

(it mm mrad) (tt mm mrad) (mm) (mm) 

monitor O10 6712 0.11 0.11 0.10 0.10 

Ap/p = 0.3% 0.51 0.10 0.10 0.02 0.01 0.03 

AB = 7%, AL = 7% 0.56 052 0.28 1.68 0.06 0.97 



Period Y (mm) 


9y (mrad) 


a (mm) 


e (n mm mrad) 


a 


Runl 0.84 


0.004 


4.17 


2.29 


-0.09 


Run2 -0.052 


-0.007 


4.08 


5.17 


0.14 


Run3b -0.024 


0.043 


3.97 


5.30 


0.25 


Run3c -0.366 


0.068 


4.22 


6.02 


0.34 



These errors are included in the beam orbit fit to the 
monitor measurements, and the magnitude of the resul- 
tant errors and their correlations are summarized in Ta- 
ble VI. 

TABLE VI: Systematic errors and correlations for the 
position and angle of the beam center at the baffle front 
surface. The X(Y) and 9x{9y) stand for horizontal 
(vertical) position and angle of the beam center, 
respectively. 



Period X (mm) Y (mm) Ox (mrad) Oy (mrad) corr(X,(9x) coyy(Y : 6y) 

Runl 0.38 OlS O056 CK29 0517 0T392 

Run2 0.27 0.62 0.064 0.32 0.752 0.398 

Run3b 0.28 0.58 0.064 0.29 0.614 0.386 

Run3c 0.35 0.58 0.072 0.28 0.697 0.417 



To estimate the systematic uncertainty on the width 
and divergence of the proton beam, the following error 
sources are considered: 

• The systematic error in the profile measurements 
by the proton beam monitors. 

• Effects of the momentum dispersion Ap/p, where a 
conservative estimate of 0.3% is assumed. 

• Uncertainties in the quadrupole magnet (FQ2, 3, 
and 4) field model: a conservative estimate of 7% 
is taken for the uncertainty in the magnetic field 
strength AB and the magnet effective length AL. 



B. Beam Direction 

The stability of the beam direction and intensity are 
also monitored by measuring the neutrino beam it- 
self [18]. Figure 12 shows the stability of the measured 
neutrino event rate (normalized by protons on target) 
measured with INGRID as a function of time. A de- 
crease of < 2% in the rate is observed, but at this time 
it has not been determined if the decrease arises from 
an actual reduction in the neutrino production or beam 
intensity dependence in the INGRID event analysis. Ta- 
ble VIII summarizes the averaged neutrino beam center 
measured in each run. The neutrino beam direction was 
kept within 1 mrad from the beam axis. 

The sources of the systematic error on the neutrino 
beam direction measurement are: systematic error on the 
INGRID measurement, the observed shift of the beam 
direction from the designed beam axis and the survey 
error (0.0024 mrad both in horizontal (x) and vertical 
(y) for SK, and 0.026 mrad in x and 0.038 mrad in y 
for ND280). In total, for Run 1 case, the systematic 
error on the neutrino beam direction relative to SK and 
ND280 direction is 0.34 mrad in x and 0.45 mrad in y. 
It corresponds to 0.44 mrad for the uncertainty of the 
off-axis angle. 



C. Horn Current 

In Run 1, the 2nd and 3rd horns were electrically con- 
nected in series and operated with one power supply. The 
1st horn was operated separately. During Runs 2 and 
3, all three horns were connected in series and operated 
with one power supply. All horns were usually operated 
at 250 kA except for Run 3b, during which the horns were 
operated at 205 kA. During the data taking periods, the 
monitored values of the horn current drifted within 2%. 
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FIG. 10: History of total accumulated protons and 
protons per pulse for the good quality beam data. The 
solid line shows the accumulated POT. The dot points 
show the number of protons per pulse. 
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FIG. 11: Time history of the measured muon profile 
center at MUMON in all run periods. A top and 

bottom figure shows the profile center in the 
horizontal(X) and the vertical(Y), respectively. A 
dashed line corresponds to 1 mrad at MUMON. Both 
directions are controlled within 1 mrad. 



This drift is most likely due to temperature dependence 
in the operation of the monitoring hardware, but varia- 
tions of the actual horn current have not been definitively 
ruled out. 



IV. THE NEUTRINO FLUX SIMULATION 

The prediction of the flux and spectrum of neutrinos at 
the T2K detectors (INGRID, ND280 and SK) is based on 
a simulation that begins with the primary proton beam 
upstream of the baffle and ends with the decay of hadrons 
or muons that produce neutrinos. The simulation and 
its associated uncertainties are driven by primary proton 
beam profile measurements, measurements of the horns' 
magnetic fields, and hadron production data, including 
NA61/SHINE measurements [11, 12]. 

FLUKA2008 [21, 22] is used to simulate the hadronic 
interactions in the target and baffle where the primary 
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FIG. 12: Neutrino events per 10 14 POT measured by 

INGRID (points) overlaid with mean value (dashed 
lines). The error bar represents the statistical error on 
the INGRID measurement. 



TABLE VIII: Neutrino beam direction measured by 
INGRID in each period. Each X and Y position 
includes the statistical error (first error term) and 
systematic error (second error term). 



Period X center [mrad] Y center [mrad] 
RUN1 0.009±0.052±0.336 -0.314±0.055±0.373 
RUN2 -0.028±0.027±0.333 0.050±0.030±0.374 
RUN3b -0.110±0.085±0.385 0.152±0.100±0.472 
RUN3c -0.001±0.026±0.331 0.232±0.029±0.378 



proton beam first interacts and produces the majority of 
the secondary pions, since FLUKA is found to have the 
best agreement with external hadron production data. 
Kinematic information for particles emitted from the tar- 
get is saved and transferred to the JNUBEAM simula- 
tion. JNUBEAM is a GEANT3 [23] Monte Carlo simu- 
lation of the baffle, target, horn magnets, helium vessel, 
decay volume, beam dump, and muon monitor. The ge- 
ometry of these components is based on the final mechan- 
ical drawings of the constructed beamline. JNUBEAM 
also includes the INGRID, ND280, and SK detectors, 
which are positioned in the simulation according to the 
latest survey results. Hadronic interactions arc modeled 
by GCALOR model [24, 25] in JNUBEAM. 

In JNUBEAM, particles are propagated through the 
horn magnetic field, and may interact with the horn ma- 
terial in the target station. Particles are propagated 
through the decay volume until they interact or decay. 
As described in Sec. IV B 2, neutrinos from particle de- 
cays are forced to point to SK or a randomly chosen 
point in the near detector plane. The neutrino kinematic 
variables and the probability based on the decay phase- 
space density and branching fraction are saved. The flux 
and energy spectrum are obtained from these simulated 
events by weighting according to the saved probabilities. 
In addition, the kinematic information of the initial pro- 
ton and full interaction chain producing the neutrino are 
saved to allow for re- weighting of the proton beam profile 
and hadron interactions. 
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FIG. 13: Flow diagram of the flux prediction. 
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FIG. 14: A two dimensional view of the geometrical 
set-up in the FLUKA simulation of the baffle and the 
target. 



The general simulation procedure is outlined in Fig. 13. 

A. Interaction of primary beam in the target 

The simulation of the interactions of the primary beam 
protons with the graphite of the baffle and the target core 
is performed using FLUKA2008. Incident protons are 
generated according to the measured proton beam spatial 
distribution and divergence. The kinetic energy is set to 
30 GeV. Figure 14 shows the two-dimensional projection 
of the simulated geometry. The baffle is described as a 
graphite block with the dimensions 29 x 40 x 171.145 cm 3 
and a 3.0 cm diameter cylindrical hole through the cen- 
ter. The target is modeled as a graphite cylinder 90 cm 
long and 2.6 cm in diameter. The volume inside the baf- 
fle hole and between the baffle and the target is filled 
with He gas. The generated particles are traced until 
they emerge from the model geometry and then informa- 
tion such as kinematic variables and hadron interaction 
history at that point is recorded. 

B. Tracking inside horns and helium vessel. 

Particles are generated in JNUBEAM according to 
the recorded information in the previous step, and are 
tracked through the horns and helium vessel. The 2 mm 
thick graphite tube and 0.3 mm thick titanium case sur- 
rounding the target core are also modeled in JNUBEAM. 
The interaction of generated particles with the materials 
in JNUBEAM is modeled by GCALOR. 

1. Horn magnetic field 

As explained in Sec. II Bl, a toroidal magnetic field 
is generated in the horns. The field strength varies as 
1/r, where r is the distance from the horn axis. Since a 



low frequency pulsed current (3.6 ms full width) is loaded 
into the horn, the skin effect is small (the estimated skin 
depth is approximately 5 mm while the thickness of the 
inner conductor is 3 mm.). Therefore, we assume that 
the current flows in the conductor uniformly. On this 
assumption, the magnetic field at radius r in the inner 
conductor is calculated with Ampere's Law as: 



where /in is the magnetic permeability, / is the current 
and a and b are, respectively, the inner and outer radii 
of the inner conductor. 



2. Neutrino production 

The particles are tracked in the helium vessel, decay 
volume, and the surrounding concrete shield including 
the beam dump until they decay or their kinetic energy 
drops below 10 MeV (at which point unstable particles 
are decayed). Decay products are also tracked except 
for neutrinos. In JNUBEAM, tt ± , K ± , and /z± de- 
cays listed in Table IX are considered as neutrino sources. 
The current best knowledge [26] on the branching ratios 
and Kf 3 (K + -> n°l + vi/K~ -> TT°l~v h I = e,(i) decay 
form factors is used. When a muon is generated from 
pion/kaon decay, its polarization information is stored. 
This polarization is then taken into account at the muon 
decays. 

In order to save computing time, when a particle de- 
cays into neutrino(s), the neutrino(s) are forced to point 
in the direction of SK or a randomly chosen point in the 
near detector planes. The neutrino energy in the center 
of mass frame is assigned based on the decay kinemat- 
ics. The neutrino is then boosted into the laboratory 
frame under the assumption that it points towards the 
desired detector, and the probability of production in the 
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selected direction is stored as an event weight. In addi- 
tion to this probability, the neutrino flavor, energy and 
ancestors' truth information are stored. The neutrino 
flux spectrum is obtained by weighting each event with 
the stored probability. For neutrinos produced with en- 
ergy less than 4 GeV, the storage of events is pre-scaled 
(and event weights are adjusted accordingly) to allow for 
sufficient statistics in the high energy tail of the flux pre- 
diction without producing prohibitively large file sets. 

TABLE IX: Neutrino-producing decay modes 
considered in JNUBEAM and their branching ratio in 
percentage. Decay modes for and D e are omitted in 

this table. The 7r~, K~ and n~ modes are charge 
conjugates of the n + , K+ and ^+ modes, respectively. 



Particle Decay Products Branching Fraction (%) 
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The hadrons are labeled as secondary hadrons if they 
are produced in interactions of the original protons, and 
tertiary hadrons if they are produced by interactions of 
hadrons other than the original proton. The breakdown 
of the predicted flux for a given flavor by the final hadron 
in the interaction chain is shown in Table X. The v e and 
v e originating from secondary or tertiary pions are from 
subsequent muon decays. A significant fraction of the 
fluxes come from tertiary pions and kaons, so it is impor- 
tant to investigate hadron interaction data at both the 
T2K beam momentum and for lower momentum hadrons. 

TABLE X: The fraction of the neutrino flux by the final 
hadron in the interaction chain after hadron interaction 
re-weighting is applied. 



Flux percentage of each(all) fiavor(s) 

Parent v e v e 
Secondary 

tt* 60.0(55.6)% 41.8(2.5)% 31.9(0.4)% 2.8(0.0)% 

K± 4.0(3.7)% 4.3(0.3)% 26.9(0.3)% 11.3(0.0)% 

K° L 0.1(0.1)% 0.9(0.1)% 7.6(0.1)% 49.0(0.1)% 

Tertiary 

tt* 34.4(31.9)% 50.0(3.0)% 20.4(0.2)% 6.6(0.0)% 

K± 1.4(1.3)% 2.6(0.2)% 10.0(0.1)% 8.8(0.0)% 

Kl 0.0(0.0)% 0.4(0.1)% 3.2(0.0)% 21.3(0.0)% 



C. The simulation of hadronic interactions 

As discussed in Sec. IV A, the hadronic interactions in 
the target are modeled with FLUKA2008. Outside of the 
target, where GEANT3 controls the simulation, interac- 
tions are modeled with GCALOR. The chain of hadronic 
interactions for each simulated event producing a neu- 
trino is saved, and re- weighting based on hadron interac- 
tion measurements is applied to the simulated events. 

The hadron interaction data used are thin target data, 
described in Sec. IV CI, that include measurements of 
inelastic cross sections and differential hadron produc- 
tion. Unlike the case of the thin target measurements, 
particles traversing the T2K target encounter a signifi- 
cant amount of material and can undergo multiple inter- 
actions. In addition particles can also interact with the 
material outside the target. A step-by-step re-weighting 
procedure is therefore applied to the hadronic interaction 
chain in each event. The weights are applied to: 

1. differential production of n^, and K\ in the 
interactions of protons on the target materials 
(Sec. IV C 2). 

2. interaction rates for p, ^ and K ± that affect the 
rate of interactions that produce hadrons, as well 
as the attenuation of hadrons that may decay to 
produce a neutrino (Sec. IV C 3). 



1. Data used for hadronic interaction re-weighting 

The pion and kaon differential production measure- 
ments used for obtaining the T2K flux predictions are 
summarized in Table XL 

TABLE XL Differential hadron production data 
relevant for the T2K neutrino flux predictions. 



Experiment Beam Mom. (GeV/c) Target Particles 

NA61/SHINE [11] [12] 31 C tt ± , K+ 

Eichten et al. [27] 24 Be, Al, ... p, n ± , K ± 

Allaby et al. [28] 19.2 Be, Al, ... p, ty ± , K ± 

BNL-E910 [29] 6.4 - 17.5 Be tt* 



To predict the neutrino flux, T2K relies primarily on 
the measurements of pion [11] and kaon [12] yields by 
the NA61/SHINE experiment at the CERN SPS. These 
data were taken with a thin (2 cm) graphite target and 
the same proton beam energy as that of T2K. The results 
are based on the data collected in 2007 during a first, lim- 
ited statistics, run with about 6.7xl0 5 registered events. 
An additional data set, taken with the target removed, 
was used to account for the contamination by particles 
produced in interactions of the proton beam occurring 
outside the target. 

Charged particles are identified by using the measure- 
ment of the specific energy loss (dE/dx) and of the time- 
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of- flight (ToF). The combined information of ToF and 
dE/dx is needed in the 1-4 GeV/c momentum range 
where different particle species have similar values for 
their specific energy loss. A calibration of the mean 
dE/dx as a function of the momentum with an accu- 
racy of 0.1% was required to limit the systematics on the 
particle identification at the level of 1%. 

Charged pion differential production cross sections 
were measured as a function of the pion laboratory mo- 
mentum in 10 intervals of the pion laboratory polar angle 
relative to the proton beam direction, covering the range 
from to 420 mrad. The considered momenta range from 
0.2 GeV/c up to 19.6 GeV/c depending on the polar an- 
gle is illustrated in Fig. 15. For momenta above about 
7.5 GeV/c a lower limit on the polar angle is set by the 
limited detector acceptance in the forward region. The 
experimental errors, dominated by the systematic uncer- 
tainties, are discussed in Sec. VAl. 

The positive kaon production measurements were per- 
formed with a coarser data binning and for a range of the 
kinematic variables which covers about 60% of the phase 
space relevant for T2K. Limitations were imposed by the 
available statistics and by the decreased sensitivity of the 
kaon identification at larger momenta as a consequence 
of the vanishing K/p and K/tt production ratios. The 
maximum kinematic range considered is between 1.6 and 
7.2 GeV/c in momentum and between 20 and 240 mrad 
for the polar angle (Fig. 15). The experimental errors on 
the K + production cross section, mainly dominated by 
the statistical uncertainties, are discussed in Sec. VA2. 

The NA61/SHINE data cover most of the relevant 
hadron production phase space for the T2K flux, as illus- 
trated in Fig. 15, which shows the simulated momentum 
and production angle of pions and kaons that are pro- 
duced in primary proton interactions and decay to con- 
tribute to the neutrino flux at SK. More than 90% of the 
pion phase space is covered, and the K + data cover 60% 
of the kaon phase space. 

The importance of the NA61/SHINE future program 
of measurements is outlined in Sec. VA5. 

The measurements of the differential kaon production 
by Eichten et al. [27] and Allaby et al. [28] cover the for- 
ward production of high energy kaons, which has not 
been measured yet by the NA61/SHINE experiment. 
These data are used to re-weight the model predictions 
in these regions. In addition, the differential proton pro- 
duction measurements in these experiments are used to 
evaluate systematic uncertainties in secondary nucleon 
production. 

The pion production data from the BNL-E910 exper- 
iment [29] is used to evaluate systematic uncertainties 
associated with tertiary pion production. 

Measurements of the inelastic cross section for proton, 
pion, and kaon beams with carbon and aluminum targets 
are used to re-weight particle interaction rates and ab- 
sorption in the simulation. A summary of these data is 
given in Table XII. The experiments typically measure 
the inelastic cross section Oi ne i which is defined as the 



total cross section minus the elastic cross section. Some 
experiments measure <J pro d, the production cross section, 
which is defined here as: 



(2) 



Here, a qe is the quasi-elastic scattering off of individual 
nuclei. The production cross section represents the rate 
of interactions where hadrons are produced in the final 
state. 



2. Hadron differential production re-weighting 

The differential production re-weighting is evaluated 
using the differential multiplicity in the momentum, p, 
of the produced particle and its angle, 9, relative to the 
incident particle: 



dn 
dp 



(6, p in , A) = 



1 



da 



a p rod{Pin,A) dp 



(6, Pin ,A). (3) 



The cross section <J pr od(Pin, A) depends on the incident 
particle momentum, p in , and target nucleus, A. 

The differential production weight that is applied to a 
given simulated interaction that produces hadrons is the 
ratio of the production in data and simulation: 



W( Pm ,A) = 



[ dp (@>Pim ^-)]data 

[%{6,p in ,A)] MC 



(4) 



For interactions of 31 GeV/c protons on carbon that 
produce ^ or K + in the phase space covered by the 
NA61/SHINE data, the construction of the ratio in Eq. 4 
is straightforward since the differential production data 
provided is already in the form in Eq. 3, at the correct 
beam momentum, and on the correct target material. 
The weights applied to differential production in FLUKA 
simulated interactions are shown in Fig. 16. 

The re-weighting of tertiary pion production from 
nucleon interactions requires extrapolations from the 
NA61/SHINE data to lower incident nucleon momen- 
tum and other target materials, since tertiary production 
can happen in interactions within the horns (aluminum) . 
Tertiary pions can also be produced in the interactions 
of secondary neutrons, in which case data for the isospin 
symmetric reaction (p+C — > n ± +X for n+C — > n T +X) 
are used to calculate weights. The same invariance is as- 
sumed for interactions on the Al nuclei, although the 
isospin invariance of the nucleus is slightly broken. 

The scaling of differential production rates to different 
incident nucleon momenta is carried out assuming Feyn- 
man scaling [42]. The Feynman variable, x F , is defined 
as: 



Xf 



PL 



VL(inax) 



(5) 



where pi, is the longitudinal momentum of the produced 



particle in the center of mass frame and p L 



(max) 



is the 
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(a) 7T+ (b) 7T- (c) K+ 



FIG. 15: The phase space of pions and kaons contributing to the predicted neutrino flux at SK, and the regions 

covered by NA61/SHINE measurements. 

TABLE XII: Inelastic and production cross-section data used to re-weight hadron absorption probabilities. 



Data 


Beam 


Target 


Beam Momentum (GeV/c] 


Measurement 


Abrams et al. [30] 


K± 


C, Cu 


1 - 3.3 


°"incl 


Allaby et al. [31] [32] 


71--, K~ 


C, Al, . 


. 20-65 


°"incl 


Allardyce et al. [33] 


7r± 


C, Al, . 


. 0.71 - 2 


°"incl 


Bellettini et al. [34] 


P 


C, Al, . 


. 19.3, 21.5 


°"incl 


Bobchenko et al. [35] 


71--, p 


C, Al, . 


. 1.75 - 9 


°"incl 


Carroll et al. [36] 




, p C, Al, . 


. 60 - 280 


^prod 


Cronin et al. [37] 


1X~ 


C, Al 


0.73 - 1.33 


°"incl 


Chen et al. [38] 


P 


C, Al, . 


. 1.53 


°"incl 


Denisov et al. [39] 


7r±, K± 


, p C, Al, . 


.6-60 


°"incl 


Longo et al. [40] 


7T+, p 


C, Al 


3 


°"incl 


NA61/SHINE [11] 


P 


C 


31 


^prod 


Vlasov et al. [41] 




C, Al 


2 - 6.7 


a incl 



maximum allowed longitudinal momentum of the pro- 
duced particle. The weights shown in Fig. 16 are con- 
verted to an xp,pT dependence and applied to tertiary 
events based on the xf and px of those events. This re- 
weighting method assumes perfect scaling, and the sys- 
tematic effect is studied in Sec. V A using data with lower 
incident particle momenta. 

The NA61/SHINE data are also extrapolated from a 
carbon target to aluminum and used to re-weight in- 
teractions in the horn material that are modeled in 
the GEANT3 (GCALOR) simulation. The A-dependent 
scaling is carried out using a parametrization proposed 
by Bonesini et al. [43] based on works by Barton et al. [44] 
and Skubic et al. [45]: 



E 



d 3 a(A 1 ) 
dp 3 



Ai 
A 



a(x F ,PT) 



E 



d 3 a(A ) 
dp 3 



where: 



a(x F ,p T ) = (a + bx F + cx 2 F )(d + ep T ). 



(6) 



(7) 



The parameters a through e are determined by fitting 
the A-dependence in the data from Eichten et al. [27] 
and Allaby et al. [28]. Examples of the fitted A depen- 
dence for a few bins are shown in Fig. 17. In this figure, 



TABLE XIII: Parameters for material scaling. 



Bonesini et al. [43] 0.74 -0.55 0.26 0.98 0.21 
Fit to tt data 0.75 -0.52 0.23 1.0 (fixed) 0.21 

Fit to K data 0.77 -0.32 0.0 1.0 (fixed) 0.25 



the ratio of the K + production from the Al target to that 
obtained from the Be target by [27] is plotted at differ- 
ent momenta for three angular bins. The accuracy and 
precision of the scaling for the individual data points is 
discussed in Sec. V A. The fitted parameter values along 
with the values reported in [43] are listed in Table XIII. 

The NA61/SHINE pion production data are scaled to 
aluminum using the parameters in Table XIII, and the 
resulting weights applied to the production in GCALOR 
are shown in Fig. 18. The weights are calculated for 
GCALOR, since the simulation of interactions in the 
horn material is done with GEANT3. 

The re-weighting of K + and K~ production in the 
phase space not covered by NA61 /SHINE is carried out 
using the Eichten et al. [27] and Allaby et al. [28] kaon 
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FIG. 16: The differential production weights from 
NA61/SHINE data for tt+ (top), n~ (middle) and K+ 
(bottom). 





Data: 17 mrad bin 




Fit: 17 mrad bin 
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FIG. 17: Examples of the material scaling exponent a 
fit for a few angular bins in the [27] K + data. 
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production data. Since these data sets only measure the 
differential production at points that cover a small mo- 
mentum and angular range, it is necessary to interpolate 
between the data points to find the weights for interme- 
diate simulated production. A bi-cubic spline interpo- 
lation is performed to each data set separately, and the 
resulting differential production cross sections are shown 
in Fig. 19. 

Since these data sets do not include points on carbon, 
the data on Be are compared to the FLUKA prediction 
for Be at the same incident particle momentum as the 
data set. The ratios of the data and FLUKA predic- 



FIG. 18: The differential production weights for 
GCALOR from A-scaled NA61 /SHINE data for tt h 
(top), 7T~ (bottom). 



tions are evaluated and the corresponding distributions 
of weights from each data set are shown in Fig. 20. 

The weights in Fig. 20 are converted to the xp ~ Pt 
basis and applied assuming x p scaling. The Eichten data 
are used in regions covered by that data set, but not cov- 
ered by the NA61/SHINE K+ data. The Allaby data are 
used in regions covered by that data set, but not covered 
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FIG. 19: The interpolated kaon production double differential cross section measurements of Eichten et al. [27] (top) 
and Allaby et al. [28] (bottom). The markers indicate the locations of the data points. 



by either NA61/SHINE K+ data or the Eichten data. 
For regions not covered by any data, no re-weighting is 
applied and the effect is studied as part of the uncer- 
tainty, as discussed in Sec. VA2. 

The Kl multiplicity is calculated from the Eichten 
and Allaby data using a simple quark parton model 
(QPM) [46, 47]. Assuming the following conditions on 
the number densities of sea and valence quarks 



c/,. 



(8a) 
(8b) 



a relation between the number of produced (K s ) 
K + , and K~ can be established: 



N{Kl) 



N(K° S ) 



N(K+) +3N(K~) 



(9) 



After calculating the K\ production according to 
Eq. 9, the K® multiplicity is re-weighted in the same 
manner as in the case of K ± . The weights are shown in 
Fig. 21. 

Although Eq. 9 is only strictly valid for proton-proton 
collisions (n = 2), the effect of proton-neutron (n = 1) 
interactions leads to only small changes in the flux pre- 
dictions that are < 1%. It is, therefore, not considered 
at this time. 



3. Hadron interaction rate re-weighting 

In addition to re- weighting the multiplicity of hadrons 
produced in interactions of nucleons, it is necessary to re- 
weight the rate at which hadrons interact. The quantity 
that is re- weighted is the production cross section defined 
in Eq. 2. Many experiments, however, measure the in- 
elastic cross section which includes the quasi-elastic com- 
ponent. To carry out re- weighting based on the <T pro( i, 
the quasi-elastic cross section must be subtracted from 
the measurements. The quasi-elastic cross section is ex- 
trapolated from hadron+nuclcon scattering data using a 
modification of the empirical dependence derived by Bel- 
lettini et al. [34]: 



a qe = 0.8(^ p + at)A 



1/3 



(10) 



Here a^ p and af n are the elastic cross sections of the 
hadron h on the proton and neutron respectively. The 
formula is modified from Bellettini et al. to include the 
average of the elastic cross section on the proton and neu- 
tron instead of the proton only. The quasi-elastic cross 
section evaluated for proton interactions on carbon in this 
manner is shown in Fig. 22. The value of a qe — 33.1 mb is 
slightly higher than the value that NA61 /SHINE derived, 
a qe = 27.9 ± 1.5 (sys) mb [11], using Glauber model cal- 
culations. As discussed in Sec. VA4, the uncertainties 
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FIG. 20: Ratios of the interpolated charged kaon double 
differential production multiplicity measurements from 
Eichten and Allaby over the FLUKA predicted yields 
from Be for 24 GeV/c and 19.2 GeV/c proton beam 
momenta, respectively. 
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FIG. 21: Ratios of the interpolated K\ double 
differential production multiplicity derived from Eichten 
and Allaby over the FLUKA predicted yields from Be 
for 24 GeV/c and 19.2 GeV/c proton beam momenta, 
respectively. 



on the weights are conservatively set to the magnitude of 
the quasi-elastic correction used to derive the production 
cross section. 

The re-weighting of the interaction rate models the 
change in the survival probability of particles as the cross 
section changes, as well as the change in rate at a given 
interaction point. The probability that a particle with 
hadron production cross section of <J vro d travels a dis- 
tance x and interacts in the infinitesimal length Ax to 
produce hadrons is: 

P(x;a prod )= <T prod pe- x ' a * r ° dP dx' (11) 



= A 



X(7p roc ipC 



-X(7 prod p 



(12) 



Here, p is the density of nuclear targets in the mate- 
rial. When the production cross section changes, a pro d — ► 
<j' prod , the weight applied to model the change is the ratio 
of the probabilities: 



W = 



P^Xj C prod) 



_ a prod e -x(v' prod -<r prod )p_ 
&prod 



(13) 
(14) 
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FIG. 22: The elastic cross sections for protons 
scattering on protons and neutrons and the derived 
quasi-elastic cross section for a carbon target. 



The first factor in Eq. 14 is the change in interaction 
probability at that point, while the second factor is the 
attenuation of the particle flux over the distance traveled. 
For a particle that decays before interacting, the weight 
applied is: 



D. Summary of the T2K flux prediction 

The T2K flux is predicted using the simulation includ- 
ing the re-weighting of hadron interaction probabilities 
outlined here. The flux is predicted for each neutrino 
flavor at the far and near detectors. Figure 25 shows the 
flux predictions for both SK and the ND280 off-axis de- 
tector broken down by the parent particle that decays to 
the neutrino. The relative fractions of each flavor in the 
SK flux prediction for - 1.5, 1.5 - 3.0 and > 3.0 GeV 
energy ranges after re-weighting is applied are listed in 
Table XIV. The v e flux, which constitutes an irreducible 
background for the study of — > v e oscillations, ac- 
counts for less than 1% of the flux below 1.5 GeV, and the 

contamination is ~ 5%. In the intermediate (1.5 — 3.0) 
GeV energy bin, the relative fraction of increases as 
the flux becomes more dominated by forward going pions 
that are not focused, which include 7r _ that decay to v^. 
The v e fraction also increases as the contribution from 
kaon decays becomes dominant. For the high energy bin 
(> 3.0 GeV), the fraction of flux decreases, since the 
contribution from the decay of focused kaons becomes 
significant. 

TABLE XIV: The fraction of the total flux by flavor in 
bins of the neutrino energy. The fractions in parentheses 
are relative to the total flux over all neutrino energies. 



The comparison of data and simulated cross sections 
in Fig. 23 shows that FLUKA is in good agreement 
with the data, while GEANT3 (GCALOR) has signifi- 
cant disagreement at low incident particle momenta. For 
both the simulated cross section and the data, where ap- 
plicable, the quasi-elastic cross sections are subtracted. 
Therefore, no weights are applied to the FLUKA sim- 
ulation of interactions in the target, but the GEANT3 
(GCALOR) production cross sections are re-weighted to 
the FLUKA value. 



4- Hadron interaction re-weighting summary 

The hadron multiplicity re-weighting described in 
Sec. IV C 2 and the hadron interaction rate re-weighting 
described in Sec. IV C 3 are applied to the simulation to 
achieve the hadron interaction re-weighted flux predic- 
tion for T2K. The effect of the weights are seen by taking 
the energy dependent ratio of the flux with and without 
the weights applied, as shown in Fig. 24. The pion multi- 
plicity re-weighting has the largest effect at low energies, 
while the kaon multiplicity re-weighting is dominant at 
high energies. 



Energy Range (GeV) 
Flavor 0-1.5 1.5-3.0 > 3.0 

V„ 0.9363(0.8570) 0.7719(0.0391) 0.8821(0.0372) 
0.0542(0.0496) 0.1729(0.0087) 0.0795(0.0034) 
v £ 0.0085(0.0078) 0.0451(0.0023) 0.0304(0.0013) 
v £ 0.0010(0.0009) 0.0100(0.0005) 0.0080(0.0003) 



V. UNCERTAINTIES ON THE FLUX 
PREDICTION 

In this section, we discuss uncertainties on the flux 
prediction. The neutrino flux uncertainties arising from 
hadron production uncertainties (Sec.V A), proton beam 
and off-axis angle uncertainties (Sec.V B), target and 
horn alignment uncertainties (Sec.V C), horn current and 
magnetic field uncertainty (Sec.V D) are considered. 

The uncertainties on the flux prediction are studied 
by varying underlying inputs to the flux simulation (the 
hadron production model, the proton beam profile, the 
horn currents, etc.) and evaluating the effect on the pre- 
dicted flux. Two approaches are used. 

Where an error source includes a number of correlated 
underlying parameters (degrees of freedom), re- weighting 
methods are used when possible. The underlying param- 
eters are varied according to their covariance, and the 
flux prediction is re-weighted with each of N sets (typi- 
cally 500 or more) of the parameter values. The effect on 
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FIG. 23: Comparisons of u pro d measurements and the values used in the simulation (solid line for FLUKA and 
dashed line for GCALOR), for incident protons (top left) and charged pions (top right), K + (bottom left) and K~ 

(bottom right). 



the flux is evaluated by constructing a covariance matrix 
from the N re-weighted versions of the flux prediction: 

k=N 

V * = A7 E ^nom 4>iWnom (™) 
fe=l 

Here the <j>nom are the nominal flux and i specifies the 
neutrino energy bin, flavor and detector at which the flux 
is evaluated. The <j> % k are the corresponding bins of the 
k th re-weighted version of the flux. Flux uncertainties 
evaluated with this method are the hadron interaction 
uncertainties and the proton beam profile uncertainties. 

The second method for evaluating uncertainties is ap- 
plied for uncertainties represented by variations of the 
flux due to changes in a single underlying parameter. For 
these uncertainties the flux is typically re-simulated for 
variations of the parameter at ±lcr. As with the previous 
method a covariance matrix is calculated: 

Vij - ^[(^ m -<)(^ om -^)+(C oro -^)(^ oro -^)]- 

(17) 



The 4> l + and are the re-simulated flux for +lcr and 
— la variations of the underlying parameter. 

The combined uncertainty on the flux prediction is 
simply represented by the sum of the covariances from 
each independent source of uncertainty described in the 
following text. Since the flux is evaluated as a covariance 
between bins in neutrino energy, neutrino flavor, and neu- 
trino detector, the covariance between the flux prediction 
at the near and far detectors is included. The covariance 
can be used directly in an extrapolation method, or to 
calculate the uncertainty on a far-to-near ratio. 

A. Hadron interaction uncertainties 

The systematic uncertainties associated with the 
hadronic interactions come from a variety of sources. One 
of them is the experimental uncertainties in the data. An- 
other is the scaling of the differential production yields to 
different incident particle momenta (see Section IV C 2). 
In addition, the systematic effects associated with the ex- 
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FIG. 24: Ratio of the hadron interaction re-weighted flux over the nominal flux for (upper left), (upper right), 
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trapolation of the differential particle yields to different 
target materials must be considered. It is also necessary 
to estimate the contribution from the regions of particle 
production phase space not covered by the data. Finally, 
the systematic uncertainties associated with the total in- 
teraction rate (production cross section) of particles in a 
given material must be evaluated. 



1. Pion production uncertainties 

The uncertainty on the pion production multiplicity 
modeling arises from a number of sources: 

1. The uncertainty on the NA61/SHINE data used to 
re-weight the pion production multiplicity 

2. The uncertainty on the incident particle momen- 
tum scaling used to apply the NA61/SHINE data 
to interactions with lower momentum incident nu- 
cleons 

3. The uncertainty from phase space that is not cov- 
ered by the NA61/SHINE data points 



The uncertainty from the NA61/SHINE pion multi- 
plicity data points is dominated by the systematic un- 
certainties, which are described in detail elsewhere [11]. 
Figure 26 shows the total errors including statistical er- 
rors for each of the NA61/SHINE p-6 bins. The total er- 
rors are typically 5 to 10% in the most important regions 
of the phase space. The dominant sources of uncertainty 
are the correction for the feed-down from strange parti- 
cle decays and particle identification. For most sources of 
uncertainty, the systematic effect is assumed to be corre- 
lated across all NA61/SHINE bins. This is a reasonable 
assumption for the feed-down error given the correlated 
model dependence of the strange particle production. For 
the particle identification error, it is assumed that bins of 
similar momenta arc more correlated, and the systematic 
errors are modeled with a ranged correlation: 

Ci ' j =1 " (Svyc") for lPi ~ Pjl - 6 GeV/c 

(18) 

=0 for \ Pi - Pj \ > 6 GeV/c 

(19) 
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FIG. 25: The flux predictions for the SK far detector and ND280 near detector broken down by the neutrino parent 
particle type. The error bars, which are too small to be seen in most energy bins, are due to the MC statistical error. 
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FIG. 26: The fractional error on the NA61/SHINE 
measurements in each of the p — 9 bins. The gap at 7T" 1 
momentum of 1.0-1.2 GeV/c is a region with no 
NA61/SHINE data points. 



Here pi and pj are the central value for the momentum in 
each bin. The functional form with a range of 6 GeV/c 
was chosen since it gives a reasonable model for the cor- 
relations and propagates the errors conservatively. 

The NA61/SHINE data are also used to re- weight pion 
production from the interactions of nucleons in the horn 
conductor aluminum after A-dependent scaling has been 
applied. For the scaled data points, additional errors of 
5% (correlated between p — 9 bins) and 5% (uncorrelated 
between p — 9 bins) are applied to account for the scaling 
uncertainty described in Section VA2. 

The error associated with scaling the NA61/SHINE 
pion multiplicity to lower incident nucleon momenta is 
studied by carrying out an alternative method of re- 
weighting tertiary events. Proton on Be data from the 
BNL-E910 experiment at beam momenta of 12.3 GeV/c 
and 17.5 GeV/c provide an alternative source for re- 
weighting interactions at lower incident momenta. The 
data are scaled from Be to C using the method out- 
lined in Section IV C 2. Since the BNL-E910 data are 
more coarsely binned in p — 9 than the NA61/SHINE, 
data at each momentum are separately fit with the em- 
pirical parametrization developed by Bonesini et at [43] 
(BMPT) for 400 and 450 GeV/c proton on Be differential 



production data. The BMPT parametrization uses the 
radial scaling variable xr and the transverse momentum 
of the produced particle pt- The xr variable, is defined 



as: 



E c 



TP cm 
max 



(20) 



where E cm is the energy of the produced particle in the 
center of mass frame, and is the maximum energy 

that the particle can have. Taylor et at [48] found that 
the invariant cross section when parametrized in xr and 
Pt does not depend on the total center of mass energy 
yfs (so called radial scaling) for y/s > 10 GeV, while 
Feynman scaling with xp (Eq. 5) only holds at higher 
\fs. The parametrization for the production of positively 
charged pions and kaons in proton collisions on nuclei is: 



E dp3 = A (^ ~ x R ) a {l + Bx R )x^ x 



[1 



a 

—Pt 



a 

2x 



T P T ]e- a/x « PT - 

R 



(21) 



The ratio of positive to negative hadron production was 
also found to be well described by simple parametriza- 
tions: 



K 71 ") = r o(l + xrY 



r(K)=r a (l-x R y 



(22) 



(23) 



The BMPT parametrization is found to work well for 
the BNL-E910 data and provides a smooth interpolation 
of the data points. Separate fits are done for the 12.3 
GeV/c and 17.5 GeV/c data to allow breaking of the xr 
scaling in the BMPT parameters. Tuning weights are 
calculated by taking the ratio of the BMPT fit to data 
over the FLUKA prediction, and plotted in the xr — px 
space, as shown in Fig. 27. If a simulated interaction 
has an incident particle momentum between the BNL- 
E910 data sets, a linear interpolation of the weights of 
the two data sets in the incident particle momentum is 
applied. Similarly a linear interpolation is applied for 
interactions with incident particle momenta between the 
17.5 GeV/c BNL-E910 and NA61/SHINE data. The al- 
ternative method described here varies from the default 
method in that it allows for a breaking of the x scaling 
and it uses data at lower incident particle momenta to 
guide the breaking of the x scaling. The uncertainty on 
the flux is estimated by applying the two methods of re- 
weighting tertiary events and taking the difference in the 
predicted flux. 

The NA61/SHINE data cover most of the phase space 
for secondary pions that contribute to the T2K neutrino 
flux. To study the effect of pion multiplicities in the un- 
covered region, the NA61/SHINE data are fitted with 
the BMPT parametrization, which is used to extrapo- 
late the data into the uncovered region. To improve the 
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FIG. 27: The weights for FLUKA p+C interactions 
derived from the BMPT fits of BNL-E910 tt+ (top) and 
7r~ (bottom) multiplicity data with 17.5 GeV/c protons. 



FIG. 28: The BMPT fits to the NA61/SHINE pion 
production data. 



TABLE XV: Fitted BMPT parameters from 
NA61/SHINE pion data. 



BMPT 






Parameter 


n+ 


TV 


A (mb/GeV^) 


188 ± 15 


90.8 ±2.7 


B 


-0.661 ±0.379 


-1.15 ±0.07 


a 


3.40 ±0.35 


1.89 ±0.13 


P 


0.303 ±0.029 


0.461 ±0.012 


a (GeV" 1 ) 


5.37±0.14 


5.19 ±0.045 


7 


0.245 ±0.018 


0.194 ±0.005 


5 


0.799 ±0.053 


0.783 ±0.017 


ro 




1.10 ±0.031 


ri 




1.95 ±0.17 



The total uncertainty on the T2K flux prediction 
due to the modeling of pion production arises from the 
sources outlined here and the magnitude of the uncer- 
tainty on the flux is summarized in Fig. 29. The domi- 
nant source of uncertainty for the and v e flux predic- 
tions near the flux peak is from the uncertainty on the 
NA61/SHINE pion multiplicity data points. 



2. Kaon production uncertainties 

Similarly to the pion case, the uncertainty on the kaon 
production multiplicity modeling comes from a number 
of sources: 



agreement between the fits, the 7r + and 7r _ are fitted sep- 
arately. Figure 28 shows that the BMPT parametrization 
is able to reasonably fit the NA61/SHINE data and the 
parameter values are listed in Table XV. The uncer- 
tainty in the FLUKA model in the uncovered region is 
estimated as the change in the flux when the production 
is re- weighted by the BMPT fits in the uncovered region. 
In addition, the uncertainty on the flux due to the un- 
certainty on the fitted BMPT parameters is included. 



1. The uncertainty on the data used to re- weight the 
kaon production multiplicity 

2. The uncertainty on the incident particle momen- 
tum scaling used to apply the data to interactions 
with lower momentum incident nucleons 

3. The uncertainty from phase space that is not cov- 
ered by the data points 

4. The uncertainty on the kaon production outside of 
the target 
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FIG. 29: Fractional flux error due to pion production as a function of neutrino energy, for each flavor and at the 

near and far detectors. 



TABLE XVI: Summary of the fractional uncertainties 
in the kaon production data. The uncertainty in the 
overall normalization is a n- The uncertainty for a given 
data bin is <7ApA6>- The uncertainty in the normalization 
for a given angular bin is OAe- 





on 




CAB 


NA61/SHINE 


2.3% 


11-24% 




Eichten et al. 


15% 


4% 


5% 


Allaby et al. 


10% 


2-5% 


10% 



The uncertainties associated with the experimental 
kaon production data are divided into three categories. 
The first is the uncertainty in the overall normalization, 
(Tat, for each data set. This uncertainty is fully correlated 
between different momentum and angular bins and, in 
the case of Eichten et al. and Allaby et al. measure- 
ments, for K + and K~ data sets. In the second category 
are the uncertainties, OApAe, which are uncorrelated be- 
tween different data bins. These are typically statistical 
uncertainties. In the final category are the uncertain- 
ties in normalization for a given angular bin, cas ■ These 
are treated as fully correlated for all momentum bins in 
each A9 for both K + and K~ data, but are taken to be 
uncorrelated between different angular bins. The mag- 
nitudes of the uncertainties in these three categories are 
summarized in Table XVI. 

In the case of the NA61/SHINE K+ data, the system- 
atic uncertainties (apart from the overall normalization) 
are treated as uncorrelated between different data bins. 
This is due to the fact that the dominant uncertainties for 
each bin are statistical. These uncertainties vary in the 
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FIG. 30: The BMPT fits to the NA61/SHINE K+ data. 



region 10 — 22% depending on the momentum bin, while 
the systematic uncertainties are around 4% for most of 
the bins. 

A coarse momentum and angular binning of the data 
had to be adopted for the NA61/ SHINE K+ data due 
to limited statistics. The sensitivity of the predicted 
neutrino flux to this choice of binning has been studied 
by modeling the shape of the K + production multiplic- 
ity within a given bin with the BMPT parameterization. 
The parameters in this parameterization have been deter- 
mined from a combined fit to the kaon production data 
of the NA61/SHINE, Eichten et al, and Allaby et al (see 
Figs. 30 - 32). The change in the predicted flux when 
such shape information is included is treated as an addi- 
tional systematic uncertainty. 
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31: The BMPT fits to the kaon data of Eichten et 
al. [27]. 
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32: The BMPT fits to the kaon data of Allaby et 
al. [281. 



Since the data of Eichten et al. and Allaby et al. used 
in the re- weighting of the kaon multiplicity are for Be tar- 
gets, uncertainties due to scaling of the differential yields 
from Be to C are applied in their case. These uncer- 
tainties are estimated based on the discrepancy between 
the measurements obtained with the Al targets by these 
two experiments and the expectations derived by scal- 
ing their yields from Be to Al following the procedure 
outlined in Section IV C 2. Two types of uncertainties 
are assigned: one, <j£ ias , based on the average discrep- 
ancy observed for all the data bins and the other, cr^ MS 
based on the RMS deviation from the mean value. To 
estimate these, the distributions of the ratios of Al to 
scaled Be yields, i?Ai/Scaiod Bcj are checked for each me- 
son type and each data set and the mean and RMS are 
extracted. An example of one such distribution for the 
K + data of Eichten et al. is shown in Fig. 33. Based 
on these distributions, 5% is assigned to both crj^ ias and 
^rms ■ T ne f° rrner is treated as a normalization type of 
an uncertainty correlated between the Eichten et al. and 
Allaby et al. data sets, while the latter is applied as an 
uncorrelated uncertainty for each data point. 

The uncertainty in data scaling for different incident 
beam energies is assigned based on the change in the pre- 
dicted neutrino fluxes when an alternative scaling vari- 
able xr is adopted in place of xp- In addition, the scaling 
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FIG. 33: Distribution of the ratios of Al to scaled Be 

yields, i?Ai/Scaicd Be, for the [27] K + data. The fit 
(gaussian function), which extracts the mean and RMS 
of the distribution, is overlaid. 



is checked with the data by re-scaling the measurements 
of Allaby et al. from the 19.2 GeV/c to 24 GeV/c inci- 
dent beam momentum and then comparing them directly 
with those of Eichten et al. The discrepancy between the 
two is then included as an additional source of systematic 
uncertainty. 
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To assign the uncertainty on the parts of the produc- 
tion phase space outside of the ones covered by the data, 
the BMPT parameterization is used. Its predictions are 
compared to those of FLUKA and the difference is as- 
signed as a systematic uncertainty. The uncertainties in 
the data used to determine the parameters in the func- 
tion are also propagated. 

For the kaon production from the interactions in the Al 
around the target, the uncertainties are evaluated based 
on the comparison of the Eichten et al. measurements to 
the GCALOR predictions. GCALOR predictions are also 
compared to the NA61/SHINE K + data after those are 
re-scaled to account for the difference in the target mate- 
rials. The discrepancy between the data and the model 
predictions is treated as the systematic uncertainty. 

The different contributions to the uncertainty associ- 
ated with the kaon production are shown in Fig. 34. 



3. Secondary nucleoli production uncertainties 

Interactions of the secondary protons (neutrons) inside 
the target contribute about 16% (5%) to the neutrino 
flux. The xp — px phase space of the contributing pro- 
tons and neutrons are shown in Fig. 35a and Fig. 35b, 
respectively. 

There are two components in the proton contribution: 
one for xp < 0.9 and the other for xp > 0.9. This is 
not the case for neutrons where only those with xp < 0.9 
contribute significantly. The high xp protons are pro- 
duced in quasi-elastic scattering or scattering where soft 
pions are produced, while the contribution from xp < 0.9 
is due to hadronic production. The evaluation of the un- 
certainty for the secondary nucleon production is, conse- 
quently, separated into two regions. 

In the region with xp < 0.9 the uncertainty for the 
secondary proton production is evaluated based on the 
discrepancy between the FLUKA model and the proton 
production measurements of Eichten et al. [27] and Al- 
laby et al. [28]. As shown in Fig. 36, the FLUKA model 
underestimates the production in the low pr region with 
xp > 0.5. The uncertainty in the neutrino flux prediction 
is calculated by weighting the FLUKA secondary proton 
and neutron production with the ratio of data to the 
FLUKA model, and is < 10%, as illustrated in Fig. 38. 
This is a conservative estimate of the uncertainty since 
no constraint on the average multiplicity of nucleons is 
applied in the rc-wcighting procedure. A re-weighting 
method that requires an average nucleon multiplicity of 
1 in the region of the phase space where N/N produc- 
tion is not kinematically allowed would lead to a smaller 
estimate of the uncertainty, and will be considered in the 
future. 

In the region of proton production with xp > 0.9, the 
incident protons undergo collisions with small momen- 
tum transfer. When studying variations in the flux due 
to changes in the secondary nucleon scattering, care is 
taken to ensure that the hard nucleon multiplicity re- 



mains unity, since no additional nucleons are produced. 
Due to the lack of relevant data, a 100% uncertainty is 
assigned on the proton production multiplicity in this re- 
gion, but the effect on the flux is still relatively small 
since these nucleons are forward-going and carry most of 
the original proton momentum. 



4- Production cross-section uncertainties 

The systematic uncertainty in the production cross sec- 
tion is conservatively taken to be represented by the mag- 
nitude of the quasi-elastic correction, cr qe , applied to the 
total inelastic cross section for a given particle and at 
given beam energy. This is based on an apparent discrep- 
ancy between the cross-section measurements for protons 
of Denisov et al. [39] and those of Bellettini et al. [34] , 
Carroll et al. [36], and NA61/ SHINE [11], which may be 
indicative of the difficulty in understanding whether ex- 
periments measure the inelastic or production cross sec- 
tions. These data are plotted in Fig. 37. 

For the measurement of Bellettini et at, the quasi- 
elastic contribution of 30.4 mb [11] has been subtracted 
from the reported value of 254 mb. In addition, the 
measurements by Denisov et al. are also shown after an 
estimate of the quasi-elastic contribution has been sub- 
tracted from the reported values. The fact that after the 
subtraction the agreement between all of the four exper- 
iments is better can be interpreted as that the magni- 
tude of the discrepancy is roughly similar to the size of 
the quasi-elastic cross section. A conservative approach 
is therefore taken by using a qe as the systematic uncer- 
tainty. 

5. Summary of the hadron production uncertainties and 
prospect from future measurements 

The uncertainty on the SK flux as a function of neu- 
trino energy due to hadronic interaction uncertainties is 
shown in Fig. 38. The uncertainties at the off-axis near 
detector are similar. At low energy, the largest sources 
of uncertainty in the flux are from the secondary nu- 
cleon production and production cross sections. At high 
energy, the flux uncertainty is instead dominated by the 
experimental errors on the kaon production. 

The results of the next set of measurements from 
NA61/SHINE will reduce the overall uncertainty on the 
neutrino flux prediction. Higher statistics thin target 
data have been collected with an upgraded detector 
configuration that increases the small angle acceptance. 
These data will have reduced uncertainties and cover 
the full phase space of interest for T2K. In particular, 
the kaon production measurement will be significantly 
improved. The pion production uncertainty is already 
well controlled by the NA61/SHINE measurement, and 
the additional data will have reduced uncertainties and 
slightly larger phase space coverage. One of the major 
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FIG. 34: Fractional flux error due to kaon production as a function of neutrino energy, for each flavor and at the 

near and far detectors. 




source of systematics, the contamination of pions from 
the decays of strange particles, will be further reduced 
by the NA61 /SHINE measurement of A and K s produc- 
tion rates. 

The ultimate precision on the flux prediction will fi- 
nally be achieved through the measurements of hadron 
emission from the same (replica) target as the one used 
by T2K. With precise replica target measurements it will 
be possible to reduce the uncertainties related to the 
hadron production via reinteractions inside the target. 



NA61/SHINE has already performed a pilot analysis us- 
ing low statistics replica target data [49] to establish the 
method for re- weighting the production of pions emitted 
from the T2K target. Fig. 39 shows the neutrino flux cal- 
culated using the re-weighting of the positively charged 
pion production based on the replica target data com- 
pared to the flux obtained with the re-weighting based 
on the NA61/SHINE thin target measurements. 
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FIG. 36: Ratio of the secondary proton measurements 
from Eichten et al. and Allaby et al. and the FLUKA 
modeling of secondary protons. Each circle is a point 
from the data sets. 



taken into account. It was found that only the system- 
atic errors for the vertical center position (V) and center 
angle (9y) of the beam have a sizable effect on the neu- 
trino flux prediction. This is because these parameters 
effectively change the off-axis angle at the far detector, 
which is displaced from the beam axis predominantly in 
the vertical direction. As an example, Fig. 40 shows the 
flux change when (Y, 9y) are changed by their error sizes. 
Therefore, only these errors are considered in the evalu- 
ation of the flux uncertainty. 

A large number of flux samples were prepared with (Y, 
6y) thrown according to correlated uncertainties listed in 
Table VI. In order to avoid re-running JNUBEAM for 
these different sets of Y and By, a special sample was 
generated with a large emittance in the Y — 9y phase 
space and then weighted to reproduce each thrown pair 
of (Y, By). 

The absolute flux normalization uncertainty arises 
from the errors on the proton beam intensity measured 
by CT5, i.e. 2% as described in Sec. II A 1. 
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FIG. 37: Production cross-section measurements for 
protons on graphite targets for momenta 20-60 GeV/c. 

The data from Denisov et al. are shown with and 
without the quasi-elastic estimate subtracted since the 
quantity that is measured is ambiguous. 



B. Proton beam and off-axis angle uncertainties 

1. Proton beam systematic uncertainties 

The proton beam is generated in the simulation ac- 
cording to the measured primary proton orbit and optics 
parameters as described in Sec. Ill A. 

To study the effects of the systematics errors in the 
proton beam measurements on the neutrino flux, those 
parameters were changed within the errors listed in Ta- 
ble VI. The correlation among different parameters was 



Neutrino beam direction (off-axis angle) systematic 
uncertainties 



The neutrino beam direction is measured by INGRID 
and the results are summarized in Table VIII. 

The neutrino flux uncertainty due to the uncertainty in 
the off-axis angle is evaluated by looking at a variation of 
the neutrino flux when the SK and ND280 detectors are 
moved by 0.44 mrad (Sec. II B 2) in JNUBEAM. To save 
computational time, the neutrino flux predictions for the 
moved detectors are calculated by using the nominal flux 
predictions and rescaling the energy and weight of each 
neutrino for the moved detector position using the stored 
parent particle information. 

Figure 41 shows the variation of the neutrino flux due 
to the off-axis angle uncertainty. The flux variations at 
the SK and the ND280 off-axis detector are similar to 
each other. 



C. Target and horn alignment uncertainties 

The systematic uncertainties associated with the tar- 
get and horn alignments, discussed in Section II C, are 
summarized in Table XVII. The effects of the tar- 
get alignment were studied by rotating the target in 
JNUBEAM by 1.3 (0.1) mrad in the horizontal (verti- 
cal) plane. This configuration results in a few percent 
change in the predicted neutrino flux, which is included 
as the systematic uncertainty. 

In the case of the horn position alignment uncertain- 
ties, the effects of horn movements along each coordinate 
axis were studied. Out of the three directions only the 
uncertainty in y results in a sizable change (at a few 
percent level) in the predicted flux. Since the dominant 
contribution to this systematic uncertainty is an overall 
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FIG. 38: Fractional flux error due to hadron production uncertainties. 



o 
a. 



° 10' 
> 

CD 



© 

CM 

E 

x 

^1 



replica target, pion multiplicity error 
thin target: 
total error 

pion multiplicity error 



4 




0.5 



1.5 



2.5 



3.5 



4.5 
, (GeV) 



FIG. 39: Re-weighted flux at the far detector based 
on the NA61/SHINE thin target and replica target 
measurements. 



uncertainty in the relative alignment between the pri- 
mary beamline and the secondary beamline, it is treated 
as fully correlated between the horns. For the case of the 
horn angular alignment uncertainties, the effects of horn 
rotations in both the horizontal and vertical plane by 
0.2 mrad were studied. Only rotations of the first horn, 
however, showed any significant effect on the predicted 
neutrino flux. 
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FIG. 40: An example of the fractional change of SK 
flux when the beam center position (Y) and center 
angle (0y) measured in Run 1 are changed by la, i.e. 
set to 1.42 mm and 0.29 mrad, respectively. 



The effects of the systematic uncertainties in the target 
and horn alignments on the predicted fluxes at ND280 
and SK are summarized in Fig. 42. For neutrinos with 
energies below 7 GeV the fractional uncertainties due to 
these sources are under 3%. 
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FIG. 41: Fractional change of the flux at ND280 
(top) and SK (bottom) corresponding to the 
systematics error of the off-axis angle. 
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FIG. 42: Fractional uncertainties due to the target and 
horn alignment in the flux for ND280 (top) and SK 
(bottom). 



TABLE XVII: Summary of the horn and target 
alignment uncertainties. 
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D. Horn current and magnetic field uncertainties 

As described in Sec. II Bl, the total uncertainty of 
the horn current measurement is 1.3% and the measured 
magnetic field strength is consistent with the expected 
one within 2%. Therefore, we adopted 2% (5 kA) as 
the total uncertainty on the absolute field strength. This 
results in 2% uncertainty at most in the neutrino flux. 

The anomalous field shown at Table III is also simu- 
lated by JNUBEAM. The effect on neutrino flux is less 
than 1% for energies up to 1 GeV, and less than 4% for 
energies greater than 1 GeV. 



E. Summary of flux uncertainties 

The total flux uncertainty as a function of neutrino 
energy, as shown in Fig. 43, is dominated by the hadron 
interaction uncertainties, with a significant contribution 
to the uncertainty around the flux peak arising from the 
off-axis angle and proton beam uncertainties. Shifts in 
the off-axis angle and proton beam tend to shift the peak 
position of the flux in energy. 

The flux correlations for each neutrino flavor and en- 
ergies from 1 — 10 GeV at the near and far detector are 
shown in Fig. 44. The correlations between the near and 
far detector are significant for the flux. It is also true 
that the and v e fluxes have significant correlations 
through the hadron interaction uncertainties, indicating 
that measurements of the flux can constrain the v e 
contamination in the beam. 

Typically, the flux prediction is used in an analysis 
where it is combined with near detector data to predict 
the flux at the far detector. The uncertainty on the ratio 
of the flux predictions at the far and near detectors is 
an estimate of how the uncertainty is propagated in an 
analysis where the flux is measured at the near detector. 
As shown in Fig. 45, the uncertainty on the far/near ratio 
for the f M flux prediction is less than 2% near the flux 
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FIG. 43: Fractional flux error including all sources of uncertainties. 
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FIG. 44: Correlations of the flux for a given flavor, energy and detector. The binning on the y-axis is identical to 

the binning on the x-axis. 



peak and less than 6% for all energies. The structure 
in the far/near ratio itself arises from the fact that the 
near detector sees a line source of neutrinos and hence 
a range of off-axis angles, while the far detector sees a 
point source and only a single off-axis angle. 



VI. FLUX PREDICTION AND T2K NEUTRINO 
DATA 

The flux prediction described here, in combination 
with the NEUT [50] neutrino interaction generator, is 
used to predict the event rates at the near and far neu- 
trino detectors. Comparisons of the predictions with the 
near detector data probe the accuracy of the flux model. 



A. The INGRID direction and rate measurements 

As described in Sec. II B 2, INGRID measures the event 
rate at each neutrino detector module and reconstructs 
the neutrino beam profile [18]. The peak of the neu- 
trino beam profile is a direct measurement of the neutrino 
beam direction. 

The flux at each module is calculated, as illustrated in 
Fig. 46, which shows how the flux prediction varies 



across the horizontal modules. The neutrino interaction 
rates at each detector module are predicted using the 
flux prediction, the NEUT interaction generator, and a 
Geant4-based detector simulation. Figure 47 shows the 
predicted and measured accumulated neutrino beam pro- 
file and Table XVIII summarizes the comparison of the 
predicted and measured beam center and rate for the 
Run 1 data taking period. For this period, the proton 
beam was aimed slightly off center of the target in the 
y-direction. Therefore an offset is expected in the IN- 
GRID profile center. The predictions agree well with the 
measurements of the neutrino interaction rate and profile 
center. 

TABLE XVIII: Summary of the predicted and 
measured INGRID beam center and rate for the Run 1 
period. The systematic uncertainty only includes the 
detector efficiency uncertainty and does not include flux 
or neutrino interaction uncertainties. 



Data Prediction 
Rate [events/POT] 1.59 xlCT 14 1.53 xlO" 

Horizontal center [mrad] 0.009±0.052(stat.)±0.336(syst.) 0.064 
Vertical center [mrad] -0.314±0.055(stat.)±0.373(syst.) -0.477 
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FIG. 46: The predicted flux at each of the horizontal 
INGRID detector modules from the center module (3) 
to the edge module (0). 



B. The ND280 inclusive f M measurement 

The rate of neutrino interactions in the off-axis ND280 
near detector is predicted using the flux prediction de- 
scribed here, the NEUT neutrino interaction generator 
(version 5.1.4), and a GEANT4 Monte Carlo simulation 
of the ND280 detector. An inclusive z/ M selection is ap- 
plied to the interactions at ND280 by searching for events 
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FIG. 47: The accumulated horizontal neutrino beam 
profile reconstructed by INGRID for the Run 1 period. 

The profile of the number of events at each detector 
module is fitted with a Gaussian function. Systematic 
errors are not shown in this plot. 



with a negatively charged track originating in the fiducial 
volume of the first fine grained detector that is tracked 
by the immediately downstream time projection cham- 
ber and identified as muon-like by dE/dx. The predicted 
muon momentum distribution for this selection is com- 
pared to the measured distribution from data collected in 
Runs 1 and 2, as shown in Fig. 48. The interactions from 
neutrinos produced in pion decays tend to produce events 
with lower muon momentum (since the neutrino energy is 
typically smaller), while neutrinos from kaon decays are 
the dominant contribution for interactions with higher 
muon momenta. The predicted and measured spectra 
show good agreement within the uncertainty of the flux 
prediction, which is ~ 10% for all muon momenta. The 
ratio of the total number of measured events relative to 
the prediction is: 



R, 



data /MC 



0.956 ± 0.0U(stat.) ± OmS(flux) (24) 



Even though there are additional neutrino interaction 
model and detector systematic error uncertainties, which 
are not quoted here, the data and our prediction show 
good agreement. 



VII. CONCLUSION 

In this paper, we have described the neutrino flux pre- 
diction in the T2K experiment. The predicted neutrino 
flux and energy spectrum are based on hadron production 
data, including NA61/SHINE measurements, the proton 
beam profile measurements in T2K, and measurements 
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FIG. 48: The predicted and measured muon momentum 
spectrum at ND280 for the inclusive selection (top) and 
the fractional flux uncertainty (not including neutrino 
interaction uncertainties nor the detector systematic 
error) and deviations of the data from the prediction on 
that sample (bottom). 



of the horn magnetic fields. The systematic uncertain- 
ties on the neutrino flux are based on uncertainties from 



these experimental measurements that are inputs to the 
flux prediction. Taking into account possible correlations 
between the systematic uncertainties for different angular 
and momentum bins in the hadron production data, we 
estimate the uncertainties on the neutrino flux including 
correlations between neutrinos of different energy and at 
different detectors. The total systematic uncertainty at 
the peak energy is approximately 15 % for both the near 
and far detector where the dominant source is the hadron 
interaction uncertainties. The uncertainty on the ratio of 
the flux predictions at the far and near detectors for 
flux is less than 2 % near the flux peak and less than 6 % 
for all energies. 

The predicted flux with simulated neutrino interac- 
tions is compared with the measurements at the near 
detectors. The measurements of the beam direction and 
event rate are consistent with the prediction. 
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